using ProgressLogging
using NBodySimulator, OrdinaryDiffEq, StaticArrays
using Plots, DataFrames, StatsPlots

function setup(t)
    T = 120.0 # K
    kb = 1.38e-23 # J/K
    ϵ = T * kb # J
    σ = 3.4e-10 # m
    ρ = 1374 # kg/m^3
    m = 39.95 * 1.6747 * 1e-27 # kg
    N = 216
    L = (m*N/ρ)^(1/3)
    R = 2.25σ
    v_dev = sqrt(kb * T / m) # m/s

    _L = L / σ
    _σ = 1.0
    _ϵ = 1.0
    _m = 1.0
    _v = v_dev / sqrt(ϵ / m)
    _R = R / σ

    bodies = generate_bodies_in_cell_nodes(N, _m, _v, _L)
    lj_parameters = LennardJonesParameters(_ϵ, _σ, _R)
    pbc = CubicPeriodicBoundaryConditions(_L)
    lj_system = PotentialNBodySystem(bodies, Dict(:lennard_jones => lj_parameters));
    simulation = NBodySimulation(lj_system, (0.0, t), pbc, _ϵ/T)

    return simulation
end

setup (generic function with 1 method)

We will consider the following integrators

config(τ, at, rt) = [
    # symplectic
    (alg=VelocityVerlet, dt=τ),
    (alg=McAte2, dt=τ),
    # (alg=CalvoSanz4, dt=τ),
    # (alg=McAte5, dt=τ),
    # (alg=Yoshida6, dt=τ),
    # (alg=KahanLi8, dt=τ),
    # DPRKN
    (alg=DPRKN6, abstol=at, rtol=rt),
    # (alg=DPRKN8, abstol=at, rtol=rt),
    (alg=DPRKN12, abstol=at, rtol=rt),
    # others
    (alg=Tsit5, abstol=at, rtol=rt),
    (alg=Vern7, abstol=at, rtol=rt),
    (alg=Vern9, abstol=at, rtol=rt)
]

config (generic function with 1 method)

We will first begin with the evolution of the energy error over time for different integrators. In this case we will fix both simulation time and timestep.

function benchmark(energyerr, rts, ts, t, configs)
    simulation = setup(t)
    prob = SecondOrderODEProblem(simulation)
    for config in configs
        alg = config.alg
        sol, rt = @timed solve(prob, alg(); progress=true, progress_name="$alg", config...)
        result = NBodySimulator.SimulationResult(sol, simulation)
        ΔE(t) = total_energy(result, t) - total_energy(result, 0)
        energyerr[alg] = [ΔE(t) for t in sol.t[2:end]]
        rts[alg] = rt
        ts[alg] = sol.t[2:end]
    end
end

ΔE = Dict()
rt = Dict()
ts = Dict()
configs = config(2.3e-4, 1e-20, 1e-20)
benchmark(ΔE, rt, ts, 40., configs)

plt = plot(xlabel="Rescaled Time", ylabel="Energy error", legend=:topleft);
for c in configs
    plot!(plt, ts[c.alg], abs.(ΔE[c.alg]), label="$(c.alg), $(rt[c.alg])s", xscale=:log10, yscale=:log10)
end
plt

Now, let us consider a fixed simulation time

t = 35.0
results = DataFrame(:integrator=>String[], :runtime=>Float64[], :τ=>Float64[], :abstol=>Float64[],
    :EnergyError=>Float64[], :timesteps=>Int[], :f_evals=>Int[]);

function benchmark(energyerr, rts, bytes, allocs, nt, nf, t, configs)
    simulation = setup(t)
    prob = SecondOrderODEProblem(simulation)
    for config in configs
        alg = config.alg
        sol, rt, b, gc, memalloc = @timed solve(prob, alg();
            save_everystep=false, progress=true, progress_name="$alg", config...)
        result = NBodySimulator.SimulationResult(sol, simulation)
        ΔE = total_energy(result, t) - total_energy(result, 0)
        energyerr[alg] = ΔE
        rts[alg] = rt
        bytes[alg] = b
        allocs[alg] = memalloc
        nt[alg] = sol.destats.naccept
        nf[alg] = sol.destats.nf + sol.destats.nf2
    end
end

benchmark (generic function with 2 methods)

and change the timesteps. For adaptive integrators, we will use tolerances instead.

τs = 10 .^range(-4, -3, length=5)
ats = 10 .^range(-20, -14, length=5)
rts = 10 .^range(-20, -14, length=5)

5-element Array{Float64,1}:
 1.0e-20
 3.162277660168379e-19
 1.0e-17
 3.1622776601683793e-16
 1.0e-14

The tolerances for adaptive algoritms were chosen such that the energy error is close to the one given by the VelocityVerlet with the chosen τ.

Now, let us run the benchmarks

@progress "Variable dt" for (τ, at, rt) in zip(τs, ats, rts)
    runtime = Dict()
    ΔE = Dict()
    nt = Dict()
    nf = Dict()
    b = Dict()
    allocs = Dict()

    GC.gc()
    benchmark(ΔE, runtime, b, allocs, nt, nf, t, config(τ, at, rt))

    for (k,v) in ΔE
        push!(results, [string(k), runtime[k], at, rt, abs(v), nt[k], nf[k]])
    end
end

Let us now examine the results

results

35 rows × 7 columns

	integrator	runtime	τ	abstol	EnergyError	timesteps	f_evals
	String	Float64	Float64	Float64	Float64	Int64	Int64
1	DPRKN12	140.487	1.0e-20	1.0e-20	0.00159857	9729	177844
2	Vern7	35.046	1.0e-20	1.0e-20	0.0688734	3998	41672
3	DPRKN6	32.8794	1.0e-20	1.0e-20	0.0363926	6035	45238
4	Tsit5	18.2855	1.0e-20	1.0e-20	0.1444	3281	21717
5	VelocityVerlet	583.327	1.0e-20	1.0e-20	0.00332889	350000	700002
6	McAte2	588.858	1.0e-20	1.0e-20	0.000154077	350000	1050002
7	Vern9	161.857	1.0e-20	1.0e-20	0.0995092	11191	191682
8	DPRKN12	141.322	3.16228e-19	3.16228e-19	0.00443486	9853	180346
9	Vern7	35.0334	3.16228e-19	3.16228e-19	0.111018	3984	41642
10	DPRKN6	31.6658	3.16228e-19	3.16228e-19	0.0225537	5940	44195
11	Tsit5	18.4532	3.16228e-19	3.16228e-19	0.141024	3264	21801
12	VelocityVerlet	326.416	3.16228e-19	3.16228e-19	0.00405704	196820	393642
13	McAte2	331.795	3.16228e-19	3.16228e-19	0.00058266	196820	590462
14	Vern9	162.485	3.16228e-19	3.16228e-19	0.104935	11173	191618
15	DPRKN12	140.164	1.0e-17	1.0e-17	0.0148921	9737	177574
16	Vern7	34.5065	1.0e-17	1.0e-17	0.0841265	3925	40712
17	DPRKN6	33.174	1.0e-17	1.0e-17	0.0455213	6069	46043
18	Tsit5	18.4812	1.0e-17	1.0e-17	0.304543	3274	21813
19	VelocityVerlet	183.214	1.0e-17	1.0e-17	0.00762108	110680	221362
20	McAte2	184.581	1.0e-17	1.0e-17	0.00368266	110680	332042
21	Vern9	162.971	1.0e-17	1.0e-17	0.0881161	11225	192226
22	DPRKN12	138.089	3.16228e-16	3.16228e-16	0.00795424	9719	177430
23	Vern7	34.945	3.16228e-16	3.16228e-16	0.112498	3959	41112
24	DPRKN6	31.8677	3.16228e-16	3.16228e-16	0.0207562	5897	43922
25	Tsit5	18.0651	3.16228e-16	3.16228e-16	0.108319	3249	21465
26	VelocityVerlet	104.651	3.16228e-16	3.16228e-16	0.00370175	62240	124482
27	McAte2	103.137	3.16228e-16	3.16228e-16	0.00384686	62240	186722
28	Vern9	160.911	3.16228e-16	3.16228e-16	0.0971133	10997	188226
29	DPRKN12	142.737	1.0e-14	1.0e-14	0.0156339	9849	179842
30	Vern7	34.7386	1.0e-14	1.0e-14	0.034158	3884	40202
31	DPRKN6	31.4633	1.0e-14	1.0e-14	0.100083	5872	43768
32	Tsit5	18.4425	1.0e-14	1.0e-14	0.314547	3259	21783
33	VelocityVerlet	58.4034	1.0e-14	1.0e-14	0.00366562	35000	70002
34	McAte2	57.548	1.0e-14	1.0e-14	0.00980374	35000	105002
35	Vern9	162.873	1.0e-14	1.0e-14	0.129173	11065	190002

The energy error as a function of runtime is given by

@df results plot(:EnergyError, :runtime, group=:integrator,
    xscale=:log10, yscale=:log10, xlabel="Energy error", ylabel="Runtime (s)")

If we consider the number of timesteps instead, we obtain

@df results plot(:EnergyError, :timesteps, group=:integrator,
    xscale=:log10, yscale=:log10, xlabel="Energy error", ylabel="Number of timesteps")

and we can also see how the muber of acceleration function calls affects the runtime

@df results plot(:f_evals, :runtime, group=:integrator,
    xscale=:log10, yscale=:log10, xlabel="Number of f calls", ylabel="Runtime (s)")

---